def empirical_dist(total, min_x, max_x):
plt.clf()
sns.histplot(total, stat='probability')
ax = plt.gca()
patches = ax.patches
frequencies = [patch.get_height() for patch in patches]
bin_edges = [patch.get_xy()[0] for patch in patches]
x, y = sns.kdeplot(np.random.choice(bin_edges, size=10000, p=frequencies),clip=(min_x, max_x)).lines[0].get_data()
probabilities=y/sum(y)
return x, probabilities
cell_dict={}
median_dict={}
min_x = 0
max_x = 9
mix_number=2
max_mix=20
EB_cell_number=20000
total_cycles=10000
nFinalCells = 100000
minCellsPerCelltype = 100
RarestCelltypeRate = 0.05
#what is the minimum number of cells per genotype given RarestCelltypeRate and minCellsPerCelltype required?
# RarestCelltypeRate(x100) : 100 = minCellsPerCelltype : X
MinCellsPerGenotype = (100*minCellsPerCelltype)/(RarestCelltypeRate*100)
MinCellsPerGenotype_Rate = MinCellsPerGenotype/nFinalCells
x, probabilities = empirical_dist(GrowthRate["D0_to_D4_growthRate"].tolist(), min_x, max_x)
x1, probabilities1 = empirical_dist(GrowthRate["D4_to_D10_growthRate"].tolist(), min_x, max_x)
k=0
while mix_number <= max_mix:
cell_dict[mix_number]={}
for i in range(1,total_cycles):
generated_data = np.random.choice(x, size=mix_number, p=probabilities)
generated_data_2 = np.random.choice(x1, size=mix_number, p=probabilities1)
if sum(generated_data)>mix_number:
k+=1
generated_cell=(EB_cell_number/mix_number)*generated_data
# if sum(generated_cell*generated_data_2)>sum(generated_cell):
# cell_dict[mix_number][k]=generated_cell*generated_data_2
cell_dict[mix_number][k]=generated_cell
simulated_exp=pd.DataFrame(cell_dict[mix_number])
plt.clf()
plot_data=(simulated_exp.div(simulated_exp.sum(axis=0), axis=1)>MinCellsPerGenotype_Rate).sum(axis=0)
testplot=sns.histplot(plot_data,discrete=True,stat='probability')
testplot.set(xlabel='Number of analysable genotypes', ylabel='Frequency')
testplot.set(xticks=range(0,15,1))
plt.title('Number of genotypes mixed: '+str(mix_number))
# plt.show()
plt.savefig(figDir+'/histogram_dual_'+str(mix_number)+'.png')
median_dict[mix_number]=[plot_data.mean()-plot_data.std(),plot_data.mean(),plot_data.mean()+plot_data.std()]
mix_number+=1
print(plot_data.std())
plt.clf()
data_samplings=pd.DataFrame(median_dict).T
data_samplings['x']=data_samplings.index
# PLot empirical distribution
sns.set_style("whitegrid")
sns.set_style("ticks")
fig,ax = plt.subplots(ncols=1,nrows=1, figsize=(7, 5),dpi=100 )
sns.set_style("white")
#plt.plot(data_samplings["x"], data_samplings[1], '-', linewidth = 3, label = "Median")
sns.lineplot(data=data_samplings,x="x",y=1,color="#FDA63A",
linewidth=3, ax = ax)
sns.lineplot(data=data_samplings,x="x",y=0,color="blue",
linewidth=1.5, ax = ax)
sns.lineplot(data=data_samplings,x="x",y=2,color="blue",
linewidth=1.5, ax = ax)
plt.fill_between(data_samplings["x"], data_samplings[0], data_samplings[2] ,
alpha=0.2, color='blue',linewidth=0.0)
#plt.xlabel('pseudotime', size=30)
ax.set_ylim([0,
round(data_samplings[2].max()+1) ])
ax.set_xlim([0, data_samplings["x"].max()])
ax.set_xticks([0, 2,5,10,20])
ax.set_xlabel('Mixed cell lines') # Y label
ax.set_ylabel('Recovered cell lines') # Y label
sns.despine(fig=ax.figure)
fig.savefig(figDir+"/MinGenotypeRate0.05.Simulation.svg", format='svg', bbox_inches='tight')
